clear;
N1 = 7; 
pi = acos(-1.0);
d=1.;
%a=.2/.51;% Seamouse
%a=.45/.93;  % diatom Fuhrmann
a=.45;
bbeta = 3.14159*a*a/d/d*2.0/sqrt(3.0);
% EPSA=1.33^2; % ROD DIELECTRIC
% EPSB=1.43^2;% MATRIX DIELECTRIC
EPSA=1.; % ROD DIELECTRIC
EPSB=13;% MATRIX DIELECTRIC
% The major symmetry points in the Brillouin zone are defined
gxg = 0.0; gyg = 0.0;
gxx = 0.0; gyx = 2.0*pi/sqrt(3.0);
gxj = 2.0*pi/3.0; gyj = 2.0*pi/sqrt(3.0);
% The itinerary around the brillouin zone is laid out.
gxi(1) = gxx; gyi(1) = gyx;
gxi(2) = gxg;
gyi(2) = gyg;
gxi(3) = gxj;
gyi(3) = gyj;
gxi(4) = gxx;
gyi(4) = gyx;    

 %  Construct the wave vectors used in the plane wave calculation
 
         L=0;
      gx1 = 2.0*pi; 
      gy1 = -2.0*pi/sqrt(3.0); 
      gx2 = 0.0; 
      gy2 = 2.0*pi*2.0/sqrt(3.0);
%      generate all the wavevectors within a hexagon of
%      diameter 2*n1.

for ii=-N1:N1;
   for j=-N1:N1;
       if ( (-(N1+1)<(ii+j)) & ((ii+j)<(N1+1))),
          L=L+1;
          gx(L) = gx1*ii + gx2*j;
          gy(L) = gy1*ii + gy2*j;
       end;
    end;
 end;
 

 % Construct the dielectric matrix
 
      vcel = sqrt(3.0)/2.d0*d*d;
   for ii =1:L;
       EPS(ii,ii) = bbeta*EPSA + (1.D0-bbeta)*EPSB;
   end;
   
    for ii=1:L;
         i1= ii + 1;
         for j=i1:L;
            x1 = gx(ii) - gx(j);
            x2 = gy(ii) - gy(j);
            x = sqrt(x1*x1+x2*x2)*a/d;
            EPS(ii,j) = 2.d0*pi*(EPSA-EPSB)*a*a/(vcel) * besselj(1,x)/x; 
% The matrix is symmetric
            EPS(j,ii) = EPS(ii,j);
        end;
    end
    
        
  
   
 % Use the itinerary constructed above to move around the Brillouin zone
     n=L;
 % Counter for the number of points in the itinerary
     nknt = 1;
for nk = 1:3;
        nk1 = nk+1;
      nmax= 20;
    for l=1:nmax;
         dkx =(l-1)*(gxi(nk1)-gxi(nk))/(nmax) + gxi(nk);
         dky =(l-1)*(gyi(nk1)-gyi(nk))/(nmax) + gyi(nk);
     for ii=1:n;
         a(ii,ii) = ((dkx+gx(ii))^2 + (dky+gy(ii))^2)/4.d0/pi/pi;
         b(ii,ii) = EPS(ii,ii);
         for j=ii+1:n;
             a(ii,j) = 0.d0;
             a(j,ii) = 0.d0;
             b(ii,j) = EPS(ii,j);
             b(j,ii) = EPS(ii,j);
         end
     end
% Find the eignevalues
         v=eig(a,b);
% Store the lowest 10 eigenvalues for plotting
         for li=1:10;
          rr(nknt,li)=sqrt(v(li));
          kxx(nknt)=nknt-1;
         end
 
      nknt = nknt + 1;
   end
  
end

% Plot lowest six eigenvalues

figure(5)
hold on;
plot(kxx,rr(:,1));
plot(kxx,rr(:,2));
plot(kxx,rr(:,3));
plot(kxx,rr(:,4));
plot(kxx,rr(:,5));
plot(kxx,rr(:,6));
plot(kxx,rr(:,7));
plot(kxx,rr(:,8));
plot(kxx,rr(:,9));
plot(kxx,rr(:,10));
ylabel('a/\lambda'); 
xlabel('wave vector itinerary');
legend('E polarization')


